1 Packages
library(dplR) # dendrochronology
library(tidyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(Hmisc) #for correlation matrix## Loading required package: lattice## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(corrplot)## corrplot 0.84 loaded2 Read in Data
Locations of Chronologies
2.1 Cook - Longfellow Trail - PIST - ITRDB PA008
NOAA Study Page: https://www.ncdc.noaa.gov/paleo-search/study/2993
Chronology file name : PA008.CRN Measurement file name : PA008.RWL
Date checked : 02DEC94 Technician’s name : MARIETTE SEKLECKI Supervisor’s name : HENRI D. GRISSINO-MAYER Beginning year : 1679 Ending year : 1981 Principal investigators: ED COOK Site name : LONGFELLOW TRAIL Site location : PENNSYLVANIA, USA Species information : PIST WHITE PINE Latitude : 4120N Longitude : 7912W Elevation : N/A M Series intercorrelation: 0.658 Avg mean sensitivity : 0.211 Avg standard deviation : 0.591 Avg autocorrelation : 0.856 Number dated series : 24 Segment length tested : 50 Number problem segments: 0 Pct problem segments : 0.00
Are there obvious misdated series? NO Number possible misdated series : N/A Percent misdated series : 0.00
Comments: NO ELEVATION GIVEN EXELLENT HIGH QUALITY CHRONOLOGY
pa008.crn <- read.crn("data/pa008.crn")## There appears to be a header in the crn file
## There is 1 seriescrn.plot(pa008.crn, add.spline = TRUE, nyrs=30)
pa008.rwl <- read.crn("data/pa008.rwl")## There does not appear to be a header in the crn file## Using cols 4-10 for decade field## There is 1 seriesrwl.report(pa008.rwl)## Number of dated series: 2
## Number of measurements: 152
## Avg series length: 19001.5
## Range: 19002
## Span: 2110 - 21111
## Mean (Std dev) series intercorrelation: -0.008042467 (0)
## Mean (Std dev) AR1: -0.3685 (0.1548564)
## -------------
## Years with absent rings listed by series
## Series 058 -- 19111 20211
## 2 absent rings (1.316%)
## -------------
## Years with internal NA values listed by series
## Warning: Using internal NA values is not standard practice and can break dplR
## Series samp.depth -- 3111 3211 5114 5115 5117 6114 6115 6214 6215 8214 8215 19114 191152.2 Guyette - Dividing Lake - PIST - ITRDB CANA127
NOAA Study Page: https://www.ncdc.noaa.gov/paleo-search/study/3448
Chronology file name : CANA127.CRN Measurement file name : CANA127.RWL Date checked : 17JUL96 Technician’s name : MARIETTE SEKLECKI Supervisor’s name : HENRI D. GRISSINO-MAYER Beginning year : 1662 Ending year : 1994 Principal investigators: R.P. GUYETTE Site name : DIVIDING LAKE Site location : ONTARIO, CANADA Species information : PIST EASTERN WHITE PINE Latitude : 4525N Longitude : 7836W Elevation : 500 M Series intercorrelation: 0.572 Avg mean sensitivity : 0.209 Avg standard deviation : 0.648 Avg autocorrelation : 0.828 Number dated series : 47 Segment length tested : 50 Number problem segments: 34 Pct problem segments : 8.15
Are there obvious misdated series? YES Number possible misdated series : 3 Percent misdated series : 6.38 Do they affect chronology quality? N/A Recommend withhold from ITRDB? N/A
Comments: #4 DLW04B, #5 DLW05B, #40 DLW19A ARE MISDATED #8 DLW09A 1736-1849, #32 LDW09B 1741-1824 HAVE LOW CORRELATIONS DATA SET COULD BE IMPROVED
cana127.crn <- read.crn("data/cana127.crn")## There appears to be a header in the crn file
## There is 1 seriescrn.plot(cana127.crn, add.spline = TRUE, nyrs=30)
cana127.rwl <- read.rwl("data/cana127.rwl")## Attempting to automatically detect format.
## Assuming a Tucson format file.
## There appears to be a header in the rwl file
## There are 47 series
## 1 dlw01a 1764 1994 0.01
## 2 dlw02a 1716 1994 0.01
## 3 dlw03a 1762 1994 0.01
## 4 DLW04B 1688 1994 0.01
## 5 dlw05b 1771 1994 0.01
## 6 dlw07a 1889 1994 0.01
## 7 dlw08a 1788 1994 0.01
## 8 DLW09A 1736 1994 0.01
## 9 dlw10a 1798 1994 0.01
## 10 dlw11a 1741 1994 0.01
## 11 dlw12a 1704 1994 0.01
## 12 dlw13a 1743 1994 0.01
## 13 dlw14a 1662 1994 0.01
## 14 dlw16a 1676 1991 0.01
## 15 dlw15a 1830 1994 0.01
## 16 dlw17a 1730 1992 0.01
## 17 dlw18b 1777 1994 0.01
## 18 dlw19b 1862 1994 0.01
## 19 dlw20b 1805 1994 0.01
## 20 dlw23b 1762 1994 0.01
## 21 dlw24b 1775 1994 0.01
## 22 dlw25b 1777 1994 0.01
## 23 dlw26b 1790 1994 0.01
## 24 dlw27b 1757 1994 0.01
## 25 dlw029 1800 1994 0.01
## 26 dlw01b 1770 1994 0.01
## 27 dlw03b 1762 1994 0.01
## 28 dlw05a 1773 1994 0.01
## 29 dlw06a 1757 1994 0.01
## 30 dlw07b 1880 1994 0.01
## 31 dlw08b 1766 1994 0.01
## 32 dlw09b 1741 1994 0.01
## 33 dlw11c 1687 1994 0.01
## 34 dlw12b 1665 1992 0.01
## 35 dlw14b 1669 1994 0.01
## 36 dlw15b 1860 1994 0.01
## 37 dlw16b 1693 1994 0.01
## 38 dlw17b 1692 1992 0.01
## 39 dlw18a 1765 1940 0.01
## 40 dlw19a 1784 1994 0.01
## 41 DLW20A 1763 1994 0.01
## 42 dlw22a 1830 1994 0.01
## 43 dlw24a 1850 1994 0.01
## 44 DLW25A 1764 1994 0.01
## 45 dlw26a 1707 1994 0.01
## 46 dlw27a 1752 1994 0.01
## 47 dlw28a 1830 1994 0.01rwl.report(cana127.rwl)## Number of dated series: 47
## Number of measurements: 10839
## Avg series length: 230.617
## Range: 333
## Span: 1662 - 1994
## Mean (Std dev) series intercorrelation: 0.5399936 (0.1049041)
## Mean (Std dev) AR1: 0.815617 (0.1225253)
## -------------
## Years with absent rings listed by series
## None
## -------------
## Years with internal NA values listed by series
## Noneseg.plot(cana127.rwl)
2.3 Guyette - Dividing Lake Aquatic - PIST - ITRDB CANA148
NOAA Study Page: https://www.ncdc.noaa.gov/paleo-search/study/3449
Chronology file name : CANA148.CRN Measurement file name : CANA148.RWL Date checked : 14MAR05 Checked by : H. ADAMS AND J. LUKAS Beginning year : 950 Ending year : 1993 Principal investigators: R.P. GUYETTE, B. COLE
Site name : DIVIDING LAKE AQUATIC Site location : CANADA
Species information : PIST EASTERN WHITE PINE Latitude : 4524 Longitude : -07836 Elevation : 400M
Series intercorrelation: 0.561 Avg mean sensitivity : 0.182 Avg standard deviation : 0.589 Avg autocorrelation : 0.862 Number dated series : 50 Segment length tested : 50
Number problem segments: 29 Pct problem segments : 8.22
cana148.crn <- read.crn("data/cana148.crn")## There appears to be a header in the crn file
## There is 1 seriescrn.plot(cana148.crn, add.spline = TRUE, nyrs=30)
cana148.rwl <- read.rwl("data/cana148.rwl")## Attempting to automatically detect format.
## Assuming a Tucson format file.
## There appears to be a header in the rwl file
## There are 50 series
## 1 DLA015 1648 1908 0.01
## 2 DLA025 1658 1974 0.01
## 3 DLA014 1675 1944 0.01
## 4 DLA020 1675 1947 0.01
## 5 DLA024 1843 1945 0.01
## 6 DLA011 1680 1920 0.01
## 7 DLA007 1551 1668 0.01
## 8 DLA013 1665 1874 0.01
## 9 DLA023 1699 1906 0.01
## 10 DLA008 1737 1915 0.01
## 11 DLA018 1665 1873 0.01
## 12 DLA012 1707 1823 0.01
## 13 DLA005 1835 1992 0.01
## 14 DLA03 1793 1983 0.01
## 15 DLA022 1709 1912 0.01
## 16 dla100 1531 1756 0.01
## 17 DLA097 1719 1967 0.01
## 18 DLA069 1568 1692 0.01
## 19 DLA109 1555 1668 0.01
## 20 DLA067 1637 1782 0.01
## 21 DL122B 1555 1780 0.01
## 22 DLA087 1684 1800 0.01
## 23 DLA044 1552 1709 0.01
## 24 DLA098 1776 1917 0.01
## 25 DLA060 1398 1574 0.01
## 26 DLA041 1625 1895 0.01
## 27 DLA047 1687 1804 0.01
## 28 DLA111 1552 1688 0.01
## 29 DLA104 1477 1725 0.01
## 30 DLA107 1675 1876 0.01
## 31 DLA051 1884 1993 0.01
## 32 DLA066 950 1152 0.01
## 33 DLA126 1650 1822 0.01
## 34 DLA101 1662 1869 0.01
## 35 DLA064 1320 1436 0.01
## 36 DLA105 1665 1810 0.01
## 37 DLA075 1705 1938 0.01
## 38 DLA128 1800 1954 0.01
## 39 DLA068 1042 1319 0.01
## 40 DLA123 1750 1948 0.01
## 41 DLA130 1605 1831 0.01
## 42 DLA057 1507 1701 0.01
## 43 DLA052 1754 1944 0.01
## 44 DLA050 1636 1803 0.01
## 45 DLA056 1568 1734 0.01
## 46 DLA131 1730 1841 0.01
## 47 DLA63 1392 1493 0.01
## 48 DLA122 1517 1650 0.01
## 49 DLA65B 1322 1442 0.01
## 50 DLA088 1388 1550 0.01rwl.report(cana148.rwl)## Number of dated series: 50
## Number of measurements: 9119
## Avg series length: 182.38
## Range: 1044
## Span: 950 - 1993
## Mean (Std dev) series intercorrelation: 0.4970208 (0.1131863)
## Mean (Std dev) AR1: 0.84128 (0.09196873)
## -------------
## Years with absent rings listed by series
## None
## -------------
## Years with internal NA values listed by series
## Noneseg.plot(cana148.rwl)
2.4 Wyse - Platte Plains - PIST - ITRDB MI015
NOAA Study Page: https://www.ncdc.noaa.gov/paleo/study/5358
Platte Plains, Michigan Pinus strobus - MI015 Additional Site Information Thomas C. Wyse, P. Charles Goebel
Purpose of Collection:
Forest disturbance history and recruitment patterns.
Cores are stored at: Forest Ecosystem Restoration & Ecology School of Natural Resources OARDC The Ohio State University 1680 Madison Avenue Wooster, OH 44691-4096
Description: The site is a sandy lake plain adjacent to Lake Michigan in northwest lower Michigan.
These cores were used with others to investigate forest disturbance history and recruitment patterns.
All ring width measurements were done with WinDENDRO. Cores were visually cross dated, and the cross dating was verified with COFECHA.
mi015.rwl <- read.rwl("data/mi015.rwl")## Attempting to automatically detect format.
## Assuming a Tucson format file.
## There appears to be a header in the rwl file
## There are 21 series
## 1 36-1 1920 2002 0.01
## 2 13-3 1912 2002 0.01
## 3 19-3 1937 2002 0.01
## 4 28-1 1933 2002 0.01
## 5 5-1 1937 2002 0.01
## 6 10-4 1928 2002 0.01
## 7 26-3 1929 2002 0.01
## 8 12-4 1962 2002 0.01
## 9 1-4 1950 2002 0.01
## 10 2-1 1936 2002 0.01
## 11 32-3 1916 2002 0.01
## 12 12-1 1958 2002 0.01
## 13 34-2 1943 2002 0.01
## 14 6-4 1987 2002 0.01
## 15 33-1 1917 2002 0.01
## 16 9-2 1916 2002 0.01
## 17 7-2 1973 2002 0.01
## 18 14-4 1973 2002 0.01
## 19 28-3 1922 2002 0.01
## 20 6-2 1928 2002 0.01
## 21 8-1 1984 2002 0.01rwl.report(mi015.rwl)## Number of dated series: 21
## Number of measurements: 1302
## Avg series length: 62
## Range: 91
## Span: 1912 - 2002
## Mean (Std dev) series intercorrelation: 0.4915569 (0.1124546)
## Mean (Std dev) AR1: 0.7149524 (0.1194401)
## -------------
## Years with absent rings listed by series
## Series 19-3 -- 2002
## 1 absent rings (0.077%)
## -------------
## Years with internal NA values listed by series
## Noneseg.plot(mi015.rwl)
plot(mi015.rwl, plot.type="spag")
Detrend all series (with ModNegExp)
mi015.rwi <- detrend(rwl = mi015.rwl, method = "ModNegExp")Build Chronology
mi015.crn <- chron(mi015.rwi)
plot(mi015.crn, add.spline=TRUE, nyrs=30)
3 Prepare Data
3.1 Extract
pa008 <- tibble::rownames_to_column(pa008.crn, "year") #add years from row names
pa008$year <- as.numeric(pa008$year) #year as a numeric
names(pa008)[names(pa008)=="STNDRD"] <- "pa008" #re-name values as data files
pa008 <- select(pa008, -samp.depth) #remove sample depth
str(pa008)## Classes 'crn' and 'data.frame': 303 obs. of 2 variables:
## $ year : num 1679 1680 1681 1682 1683 ...
## $ pa008: num 0.755 0.939 1.218 0.975 1.224 ...cana127 <- tibble::rownames_to_column(cana127.crn, "year") #add years from row names
cana127$year <- as.numeric(cana127$year)
names(cana127)[names(cana127)=="dlwSTD"] <- "cana127" #re-name values as data files
cana127 <- select(cana127, -samp.depth)
str(cana127)## Classes 'crn' and 'data.frame': 333 obs. of 2 variables:
## $ year : num 1662 1663 1664 1665 1666 ...
## $ cana127: num 0.949 1.161 1.788 1.552 0.974 ...cana148 <- tibble::rownames_to_column(cana148.crn, "year") #add years from row names
cana148$year <- as.numeric(cana148$year)
names(cana148)[names(cana148)=="DLAQUS"] <- "cana148" #re-name values as data files
cana148 <- select(cana148, -samp.depth)
str(cana148)## Classes 'crn' and 'data.frame': 1044 obs. of 2 variables:
## $ year : num 950 951 952 953 954 955 956 957 958 959 ...
## $ cana148: num 0.872 0.961 0.926 0.621 0.548 ...mi015 <- tibble::rownames_to_column(mi015.crn, "year") #add years from row names
mi015$year <- as.numeric(mi015$year)
names(mi015)[names(mi015)=="xxxstd"] <- "mi015" #re-name values as data files
mi015 <- select(mi015, -samp.depth)
str(mi015)## Classes 'crn' and 'data.frame': 91 obs. of 2 variables:
## $ year : num 1912 1913 1914 1915 1916 ...
## $ mi015: num 0.492 0.548 0.477 0.388 0.497 ...3.2 Filter for Common Date
pa008 <- filter(pa008, year > 1678)
cana127 <- filter(cana127, year > 1678)
cana148 <- filter(cana148, year > 1678)
mi015 <- filter(mi015, year > 1678)3.3 Join Data
all.crn <- pa008 %>% #in matrix format
full_join(cana127, by = "year") %>%
full_join(cana148, by = "year") %>%
full_join(mi015, by = "year")
tidy.crn <- gather(all.crn, "source", "std", 2:5) #in tidy format4 Analyze Data
4.1 Plot
Since 1850
ggplot(tidy.crn, aes(x=year, y=std, col = source)) +
geom_line(alpha = 0.5) +
#facet_grid(source ~ .) +
scale_x_continuous(limits = c(1850, 2010)) +
theme_classic() ## Warning: Removed 784 rows containing missing values (geom_path).
Since 1968
ggplot(tidy.crn, aes(x=year, y=std, col = source)) +
geom_line(alpha = 0.5) +
theme_classic() ## Warning: Removed 271 rows containing missing values (geom_path).
Without mi015
ggplot(data=(tidy.crn %>% filter(source != "mi015")), aes(x=year, y=std, col = source)) +
geom_line(alpha = 0.5) +
theme_classic() ## Warning: Removed 38 rows containing missing values (geom_path).
4.2 Correlation Matrix
Using description & code from: http://www.sthda.com/english/wiki/correlation-matrix-a-quick-start-guide-to-analyze-format-and-visualize-a-correlation-matrix-using-r-software
4.2.1 All Chronologies (1912-1993)
matrix <- filter(all.crn, year > 1911, year < 1994)
matrix <- select(all.crn, -year)cor.crn <- rcorr(as.matrix(matrix))
cor.crn## pa008 cana127 cana148 mi015
## pa008 1.00 0.26 0.23 -0.23
## cana127 0.26 1.00 0.71 0.22
## cana148 0.23 0.71 1.00 0.03
## mi015 -0.23 0.22 0.03 1.00
##
## n
## pa008 cana127 cana148 mi015
## pa008 303 303 303 70
## cana127 303 316 315 83
## cana148 303 315 315 82
## mi015 70 83 82 91
##
## P
## pa008 cana127 cana148 mi015
## pa008 0.0000 0.0000 0.0557
## cana127 0.0000 0.0000 0.0496
## cana148 0.0000 0.0000 0.8008
## mi015 0.0557 0.0496 0.8008# flattenCorrMatrix Function
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
row = rownames(cormat)[row(cormat)[ut]],
column = rownames(cormat)[col(cormat)[ut]],
cor =(cormat)[ut],
p = pmat[ut]
)
}Results
as_tibble(flattenCorrMatrix(cor.crn$r %>% round(2), cor.crn$P %>% round(3)))## # A tibble: 6 x 4
## row column cor p
## <fct> <fct> <dbl> <dbl>
## 1 pa008 cana127 0.26 0
## 2 pa008 cana148 0.23 0
## 3 cana127 cana148 0.71 0
## 4 pa008 mi015 -0.23 0.056
## 5 cana127 mi015 0.22 0.05
## 6 cana148 mi015 0.03 0.801Correlation plot, blank if correlation is insignificant
corrplot(cor.crn$r, type="upper", order="hclust",
p.mat = cor.crn$P, sig.level = 0.05, insig = "blank", tl.col = "black")
4.2.2 Older Chronologies (-mi015, 1679-1982)
matrix2 <- filter(all.crn, year < 1982)
matrix2 <- select(all.crn, -year, -mi015)cor.crn2 <- rcorr(as.matrix(matrix2))
cor.crn2## pa008 cana127 cana148
## pa008 1.00 0.26 0.23
## cana127 0.26 1.00 0.71
## cana148 0.23 0.71 1.00
##
## n
## pa008 cana127 cana148
## pa008 303 303 303
## cana127 303 316 315
## cana148 303 315 315
##
## P
## pa008 cana127 cana148
## pa008 0 0
## cana127 0 0
## cana148 0 0Results
as_tibble(flattenCorrMatrix(cor.crn2$r %>% round(2), cor.crn2$P))## # A tibble: 3 x 4
## row column cor p
## <fct> <fct> <dbl> <dbl>
## 1 pa008 cana127 0.26 0.00000421
## 2 pa008 cana148 0.23 0.0000458
## 3 cana127 cana148 0.71 0Correlation plot, blank if correlation is insignificant
corrplot(cor.crn2$r, type="upper", order="hclust",
p.mat = cor.crn2$P, sig.level = 0.05, insig = "blank", tl.col = "black")